#This script will intersect any number of ASCII layers stored in a single directory with XY data
rm(list=ls())
in.dir="C:/R/In"
ASCII.dir="C:/R/In/ASCII"
out.dir="C:/R/Out"
setwd(in.dir)
library("SDMTools")


#First import XY data in UTM (east/north).  Make sure these files only contain XY data or it will fuck up the extract.data command later

raw.sample=read.csv("outsample.intersect.csv", header=T)
sample.out=xy.sample=data.frame(raw.sample$raw.sample.east, raw.sample$raw.sample.north)
raw.georef=read.csv("georefxy.csv", header=T)
georef.out=xy.georef=data.frame(raw.georef$east, raw.georef$north)

#List all ASCII files in ASCII.dir

ASCII.names=list.files(ASCII.dir)
#ASCII.names=c("slope.asc","latdec.asc")
num=2

#Create a loop that reads in one ASCII file at a time and intersects it with XY sample data

for(ASCII in ASCII.names)

  {
  
  num=num+1
  t.asc = read.asc(paste(ASCII.dir,ASCII, sep="/"))
  t.int = extract.data(xy.sample,t.asc)
  #t.out = data.frame(xy.sample,t.int)
  sample.out=cbind(sample.out, t.int)
  ncol(sample.out)
  names(sample.out)[num]=ASCII
  rm(t.asc)
  
  }
  
#Write out data to a .csv file
data.out=data.frame(raw.data, sample.out)  
write.csv(x=data.out,file=paste(out.dir, "sample.intersect.csv", sep=""), row.names=F)
nrow(sample.out)

#Create a loop that reads in one ASCII file at a time and intersects it with XY georef data

for(ASCII in ASCII.names)

  {

  num=num+1
  t.asc = read.asc(paste(ASCII.dir,ASCII, sep="/"))
  t.int = extract.data(xy.georef,t.asc)
  t.out = data.frame(xy.georef,t.int)
  georef.out=cbind(georef.out, t.out[,3])
  ncol(georef.out)
  names(georef.out)[num]=ASCII
  rm(t.asc)

  }
  
# Write out data as a .csv

write.csv(x=georef.out,file=paste(out.dir, "georef.intersect.csv", sep=""), row.names=F)

#Insert final command so that the HPC runs the last write.csv

ncol(georef.out)
  

  



